###Diplomatic Statements and the Strategic Use of Terrorism in Civil Wars###
###Levy, Dudley, Chen, Siegel###
###Journal of Conflict Resolution###
###Functions###

## functions for linear model
LM_coef <- function(modResults, data, vars = NULL) {
  modelcoef=lapply(modResults, 
                   function(x) FUN=coef(x))
  modelse = lapply(modResults, 
                   function(x) FUN=sqrt(diag(vcov(x))))
  
  varnames =lapply(modResults, function(x) FUN = names(coef(x)))
  
  dfplot <- list()
  
  for (i in 1:length(modelcoef)){
    ad = 0.5 / modelse[[i]]
    modelcoef[[i]] = modelcoef[[i]]* ad
    modelse[[i]] = modelse[[i]]*ad 
    
    ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
    yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
    ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
    yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
    dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                             ylo90, yhi90,
                             Model = paste('Model',seq_along(modelcoef)[i]))
  }
  coefficientdata = do.call(rbind,dfplot)
  colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                                 "Low95CI", "High95CI", 
                                 "Low90CI", "High90CI",
                                 "Model.Name")
  
  df <- coefficientdata %>%
    dplyr::mutate(Variables = as.character(Variables)) %>%
    dplyr::arrange(desc(Variables))
  
  if(!is.null(vars)){
    df <-  df[grep(paste(vars, collapse = "|"), df$Variables),]
    
    #re-order the variables: if factor, the order of factor levels is used; if character, an alphabetical order ist used
    df$Variables <- factor(df$Variables, levels = vars)
  }
  
  
  p = ggplot(df, aes(colour =Model.Name)) + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_linerange(aes(x = Variables, ymin = Low90CI,
                       ymax = High90CI),
                   lwd = 1, position = position_dodge(width = 1/2)) +
    geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                        ymax = High95CI, shape = Model.Name),
                    lwd = 0.5, position = position_dodge(width = 1/2),
                    fill = "WHITE") +
    coord_flip() + 
    theme_bw() + xlab("") + ylab("Rescaled Coefficient") + 
    theme(legend.position="bottom",
          legend.title=element_blank(),
          axis.text = element_text(size=14),
          text = element_text(size=14),
          plot.title = element_text(hjust = .5, size = 14, face = "bold"))
  return(p)   
}  


## functions
MLM_coef <- function(modResults, data, vars = NULL) {
  modelcoef=lapply(modResults, 
                   function(x) FUN=fixef(x))
  modelse = lapply(modResults, 
                   function(x) FUN=sqrt(diag(vcov(x))))
  
  varnames =lapply(modResults, function(x) FUN = names(fixef(x)))
  
  dfplot <- list()
  
  for (i in 1:length(modelcoef)){
    ad = 0.5 / modelse[[i]]
    modelcoef[[i]] = modelcoef[[i]]* ad
    modelse[[i]] = modelse[[i]]*ad 
    
    ylo95 =modelcoef[[i]] - qt(.975, nrow(data))*(modelse[[i]])
    yhi95 =modelcoef[[i]] + qt(.975, nrow(data))*(modelse[[i]])
    ylo90 =modelcoef[[i]] - qt(.95, nrow(data))*(modelse[[i]])
    yhi90 =modelcoef[[i]] + qt(.95, nrow(data))*(modelse[[i]])
    dfplot[[i]] = data.frame(varnames[[i]], modelcoef[[i]], modelse[[i]], ylo95, yhi95,
                             ylo90, yhi90,
                             Model = paste('Model',seq_along(modelcoef)[i]))
  }
  coefficientdata = do.call(rbind,dfplot)
  colnames(coefficientdata) <- c("Variables", "Coefficients", "S.E",
                                 "Low95CI", "High95CI", 
                                 "Low90CI", "High90CI",
                                 "Model.Name")
  
  df <- coefficientdata %>%
    dplyr::mutate(Variables = as.character(Variables)) %>%
    dplyr::arrange(desc(Variables))
  
  if(!is.null(vars)){
    df <-  df[grep(paste(vars, collapse = "|"), df$Variables),]
    
    #re-order the variables: if factor, the order of factor levels is used; if character, an alphabetical order ist used
    df$Variables <- factor(df$Variables, levels = vars)
  }
  
  
  p = ggplot(df, aes(colour =Model.Name)) + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_linerange(aes(x = Variables, ymin = Low90CI,
                       ymax = High90CI),
                   lwd = 1, position = position_dodge(width = 1/2)) +
    geom_pointrange(aes(x = Variables, y = Coefficients, ymin = Low95CI,
                        ymax = High95CI, shape = Model.Name),
                    lwd = 0.5, position = position_dodge(width = 1/2),
                    fill = "WHITE") +
    coord_flip() + 
    theme_bw() + xlab("") + ylab("Rescaled Coefficient") + 
    theme(legend.position="bottom",
          legend.title=element_blank(),
          axis.text = element_text(size=14),
          text = element_text(size=14),
          plot.title = element_text(hjust = .5, size = 14, face = "bold"))
  return(p)   
}  



### a function to do the simulation-via observed-value for marginal effects

MLM_maginal_plot <- function(Fits, n_sim = 1000, var, val1, val2, data, modname){
  # load R pacakges      
  library(lme4)
  library(arm); library(purrr)
  library(ggplot2); library(ggthemes); library(viridis)
  library(ggridges); library(purrr); library(dplyr)
  #empty list to store marginal results       
  marginal <- list()     
  
  for (i in seq_along(Fits)){
    
    #extrat the model data
    mydata2 <- as.data.frame(model.matrix(Fits[[i]]))
    
    ### simulation based on observed values
    set.seed(321)
    draw <- sim(Fits[[i]], n.sims=1000)@fixef
    #copy dataset for marginal effect
    data0 <- data1 <- mydata2
    value <- as.numeric(quantile(mydata2[, var[i]], c(val1, val2)))
    #set the scenario for marginal effect
    data0[, var[i]] <- value[1]
    data1[, var[i]] <- value[2]
    
    #empty containers
    X1 <- as.matrix(data0)
    X2 <- as.matrix(data1)
    #store marginal effect in the list
    marginal[[modname[i]]] <- apply(apply(X2, 1, function (x) draw %*% x) - 
                                      apply(X1, 1, function (x) draw %*% x), 1, mean) #1 indicates row, 2= columns
  }
  #convert the list to a data frame for plotting
  df_plot <- map(marginal, data.frame) %>%
    map2_df(., names(.), ~mutate(.x, type = .y))
  #rename first variable
  names(df_plot)[1] <- "value"
  
  ##change the order of the variable
  df_plot$type <- factor(df_plot$type, levels = rev(modname))
  
  #### Plotting
  fig <- ggplot(df_plot, aes(x = value, y = type, height=..density.., fill = type)) +
    geom_density_ridges(col = "grey70", scale = 1.2, show.legend = F) +
    scale_fill_viridis(discrete = TRUE) +
    geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
    theme_ridges(font_size = 13, grid = F, center_axis_labels = T) +
    theme(axis.title.y = element_blank(), 
          text = element_text(size=13),
          plot.title = element_text(hjust = .5, size = 15, face = "bold")) 
  ## make a return obj
  p <- list(fig, df_plot)
  #rename list
  names(p) <- c("plot", " df_plot")
  return(p)
  
}


MLM_prediction_plot <- function(Fits, n_sim = 1000, var, xlabs, val1, val2, intervals, data, modname){
  # load R pacakges      
  library(lme4)
  library(arm); library(purrr)
  library(ggplot2); library(ggthemes); library(viridis)
  library(ggridges); library(purrr); library(dplyr)
  #empty list to store marginal results       
  marginal <- list()     
  
  for (i in seq_along(Fits)){
    
    #extrat the model data
    mydata2 <- as.data.frame(model.matrix(Fits[[i]]))
    
    ### simulation based on observed values
    set.seed(321)
    draw <- sim(Fits[[i]], n.sims=1000)@fixef
    value <- seq(val1, val2, by =intervals)
    df_prediction <- array(NA, c(length(value), 4))
    #set the scenario for marginal effect
    for (j in 1:length(value)){
      data0 <- mydata2
      data0[, var[i]] <- value[j]
      #empty containers
      X1 <- as.matrix(data0)
      prediction <-  apply(apply(X1, 1, function (x) draw %*% x) , 1, mean) #1 indicates row, 2= columns
      df_prediction[j, 1] <- value[j]
      df_prediction[j, 2] <- mean(prediction)
      df_prediction[j, 3:4] <- quantile(prediction, probs = c(0.025, 0.975))
      colnames(df_prediction) <-  c("X", "mean", "lo", "hi")
    }
    #store marginal effect in the list
    marginal[[modname[i]]] <- df_prediction
  }
  #convert the list to a data frame for plotting
  df_plot <- map(marginal, data.frame) %>%
    map2_df(., names(.), ~mutate(.x, type = .y))
  #rename first variable
  names(df_plot)[1] <- "value"
  
  ##change the order of the variable
  df_plot$type <- factor(df_plot$type, levels = rev(modname))
  
  #### Plotting
  fig <- ggplot(df_plot, aes(x = value, y=mean)) +
    theme_bw() +
    geom_ribbon(aes(ymin = lo, ymax = hi), alpha=.2, color = NA) +
    geom_line(aes(y = mean, x = value), size = 0.5, alpha = 0.75) + 
    facet_wrap(~type, scales = "free") +
    scale_x_continuous(name = xlabs, breaks = seq(val1, val2, by = intervals)) +
    theme(text = element_text(size=13),
          plot.title = element_text(hjust = .5, size = 15, face = "bold"),
          strip.text.x = element_text(size = 12, color='white',
                                      angle=0),
          strip.background = element_rect(fill = "#525252", color='#525252')) 
  ## make a return obj
  p <- list(fig, df_plot)
  #rename list
  names(p) <- c("plot", " df_plot")
  return(p)
  
}



### function for interaction
MLM_interEffect_binary <- function(ModelResults, bvarname1, cvarname2,
                                   cval1, cval2, labels = c("No Pro-Rebel Resolution",
                                                            "Pro-Rebel Resolution")){
  
  ### simulation based on observed values
  set.seed(321)
  draw <- sim(ModelResults, n.sims=1000)@fixef
  
  ##set simulation
  df <- array(NA, c(2, 4))
  
  X1 <- model.matrix(ModelResults)
  X2 <- model.matrix(ModelResults)
  X1[, bvarname1] = 0
  X1[, cvarname2] = cval1
  X1[, paste(bvarname1, cvarname2, sep = ":")] = 0*cval1
  
  X2[,bvarname1] = 0
  X2[,cvarname2] = cval2
  X2[, paste(bvarname1, cvarname2, sep = ":")] = 0*cval2
  
  fd_no <- apply(apply(X2, 1, function (x) draw %*% x) - 
                   apply(X1, 1, function (x) draw %*% x), 1, mean) #1 indicates row, 2= columns
  
  X3 <- model.matrix(ModelResults)
  X4 <- model.matrix(ModelResults)
  X3[, bvarname1] = 1
  X3[, cvarname2] = cval1
  X3[, paste(bvarname1,cvarname2, sep = ":")] = 1*cval1
  
  X4[,bvarname1] = 1
  X4[,cvarname2] = cval2
  X4[, paste(bvarname1, cvarname2, sep = ":")] = 1*cval2
  
  fd_yes <- apply(apply(X4, 1, function (x) draw %*% x) - 
                    apply(X3, 1, function (x) draw %*% x), 1, mean) #1 indicates row, 2= columns
  
  df[1, 1] <- 0
  df[1, 2] <- mean(fd_no)
  df[1, 3:4] <- quantile(fd_no, probs = c(.05,.95))
  df[2, 1] <- 1
  df[2, 2] <- mean(fd_yes)
  df[2, 3:4] <- quantile(fd_yes, probs = c(.05,.95))
  
  df_plot <- as.data.frame(df)
  
  
  colnames(df_plot) <- c("type", "mean", "lo", "hi")     
  
  p <- ggplot(df_plot, aes(x=type, y=mean)) +
    geom_hline(aes(yintercept=0), linetype=2, color = "black") + 
    geom_point(size=4) + 
    geom_linerange(aes(ymin=lo, ymax=hi),alpha = 1, size = 1) + 
    theme_bw() + scale_x_discrete(limits = c(0,1), labels = labels)+
    xlab("") + ylab("Marginal Effects")+
    theme(legend.position="none",
          legend.title=element_blank(),
          axis.text = element_text(size=14),
          text = element_text(size=14),
          plot.title = element_text(hjust = .5, size = 16, face = "bold"))
  
  return(p)
}

MLM_interaction_plot <- function(ModelResults, n_sim = 1000, bvarname1, cvarname2,
                                 xlabs, val1, val2, intervals,labels = c("No Pro-Rebel Resolution",
                                                                         "Pro-Rebel Resolution")){
  # load R pacakges      
  library(lme4)
  library(arm); library(purrr)
  library(ggplot2); library(ggthemes); library(viridis)
  library(ggridges); library(purrr); library(dplyr)
  
  
  ### simulation based on observed values
  set.seed(321)
  draw <- sim(ModelResults, n.sims=1000)@fixef
  
  ##set simulation
  bvarname1_val = c(0,1)
  vcvarname2_val = seq(val1, val2, by =intervals)
  
  df_plot <- list()
  df <- array(NA, c(length(vcvarname2_val),5))
  
  for (j in 1:length(bvarname1_val)){
    for (i in 1:length(vcvarname2_val)){
      
      X1 <- model.matrix(ModelResults)
      
      X1[, bvarname1] = bvarname1_val[j]
      X1[, cvarname2] = vcvarname2_val[i]
      X1[, paste(bvarname1, cvarname2, sep = ":")] = bvarname1_val[j]*vcvarname2_val[i]
      
      fd = apply(apply(X1, 1, function (x) draw %*% x) , 1, mean)
      df[i, 1] <- bvarname1_val[j]
      df[i, 2] <- vcvarname2_val[i]
      df[i, 3] <- mean(fd)
      df[i, 4:5] <- quantile(fd, probs = c(.025,.975))
      df_plot[[j]] <- list(df)
    }
  }
  
  
  library(purrr)
  library(RColorBrewer)
  df_plot <- map(df_plot, data.frame) %>%
    map_df(., rbind)
  
  colnames(df_plot) <- c("predx", "modx",  "mean", "lo", "hi")
  ##change the order of the variable
  df_plot$predx2 <- factor(df_plot$predx, labels = labels)
  p <- ggplot(df_plot, aes(x = modx, y=mean)) +
    theme_bw() +
    geom_ribbon(aes(ymin = lo, ymax = hi), alpha=.2, color = NA) +
    geom_line(aes(y = mean, x = modx,), size = 0.5, alpha = 0.75) +
    facet_wrap(~predx2, scales = "free") +
    scale_x_continuous(name = xlabs, breaks = seq(val1, val2, by = intervals)) +
    theme(text = element_text(size=13),
          plot.title = element_text(hjust = .5, size = 15, face = "bold"),
          strip.text.x = element_text(size = 12, color='white',
                                      angle=0),
          strip.background = element_rect(fill = "#525252", color='#525252')) 
  
  return(p)
}

### a function to do the simulation-via observed-value for marginal effects

LM_maginal_plot <- function(Fits, n_sim = 1000, var, val1, val2, data, modname){
  # load R pacakges      
  library(lme4)
  library(arm); library(purrr)
  library(ggplot2); library(ggthemes); library(viridis)
  library(ggridges); library(purrr); library(dplyr)
  #empty list to store marginal results       
  marginal <- list()     
  
  for (i in seq_along(Fits)){
    
    #extrat the model data
    mydata2 <- as.data.frame(model.matrix(Fits[[i]]))
    
    ### simulation based on observed values
    set.seed(321)
    draw <- sim(Fits[[i]], n.sims=1000)
    #copy dataset for marginal effect
    data0 <- data1 <- mydata2
    value <- as.numeric(quantile(mydata2[, var[i]], c(val1, val2)))
    #set the scenario for marginal effect
    data0[, var[i]] <- value[1]
    data1[, var[i]] <- value[2]
    
    #empty containers
    X1 <- as.matrix(data0)
    X2 <- as.matrix(data1)
    #store marginal effect in the list
    marginal[[modname[i]]] <- apply(apply(X2, 1, function (x) coef(draw) %*% x) - 
                                      apply(X1, 1, function (x) coef(draw) %*% x), 1, mean) #1 indicates row, 2= columns
  }
  #convert the list to a data frame for plotting
  df_plot <- map(marginal, data.frame) %>%
    map2_df(., names(.), ~mutate(.x, type = .y))
  #rename first variable
  names(df_plot)[1] <- "value"
  
  ##change the order of the variable
  df_plot$type <- factor(df_plot$type, levels = rev(modname))
  
  #### Plotting
  fig <- ggplot(df_plot, aes(x = value, y = type, height=..density.., fill = type)) +
    geom_density_ridges(col = "grey70", scale = 1.2, show.legend = F) +
    scale_fill_viridis(discrete = TRUE) +
    geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
    theme_ridges(font_size = 13, grid = F, center_axis_labels = T) +
    theme(axis.title.y = element_blank(), 
          text = element_text(size=13),
          plot.title = element_text(hjust = .5, size = 15, face = "bold")) 
  ## make a return obj
  p <- list(fig, df_plot)
  #rename list
  names(p) <- c("plot", " df_plot")
  return(p)
  
}


LM_prediction_plot <- function(Fits, n_sim = 1000, var, xlabs, val1, val2, intervals, data, modname){
  # load R pacakges      
  library(lme4)
  library(arm); library(purrr)
  library(ggplot2); library(ggthemes); library(viridis)
  library(ggridges); library(purrr); library(dplyr)
  #empty list to store marginal results       
  marginal <- list()     
  
  for (i in seq_along(Fits)){
    
    #extrat the model data
    mydata2 <- as.data.frame(model.matrix(Fits[[i]]))
    
    ### simulation based on observed values
    set.seed(321)
    draw <- sim(Fits[[i]], n.sims=1000)
    value <- seq(val1, val2, by =intervals)
    df_prediction <- array(NA, c(length(value), 4))
    #set the scenario for marginal effect
    for (j in 1:length(value)){
      data0 <- mydata2
      data0[, var[i]] <- value[j]
      #empty containers
      X1 <- as.matrix(data0)
      prediction <-  apply(apply(X1, 1, function (x)  coef(draw) %*% x) , 1, mean) #1 indicates row, 2= columns
      df_prediction[j, 1] <- value[j]
      df_prediction[j, 2] <- mean(prediction)
      df_prediction[j, 3:4] <- quantile(prediction, probs = c(0.025, 0.975))
      colnames(df_prediction) <-  c("X", "mean", "lo", "hi")
    }
    #store marginal effect in the list
    marginal[[modname[i]]] <- df_prediction
  }
  #convert the list to a data frame for plotting
  df_plot <- map(marginal, data.frame) %>%
    map2_df(., names(.), ~mutate(.x, type = .y))
  #rename first variable
  names(df_plot)[1] <- "value"
  
  ##change the order of the variable
  df_plot$type <- factor(df_plot$type, levels = rev(modname))
  
  #### Plotting
  fig <- ggplot(df_plot, aes(x = value, y=mean)) +
    theme_bw() +
    geom_ribbon(aes(ymin = lo, ymax = hi), alpha=.2, color = NA) +
    geom_line(aes(y = mean, x = value), size = 0.5, alpha = 0.75) + 
    facet_wrap(~type, scales = "free") +
    scale_x_continuous(name = xlabs, breaks = seq(val1, val2, by = intervals)) +
    theme(text = element_text(size=13),
          plot.title = element_text(hjust = .5, size = 15, face = "bold"),
          strip.text.x = element_text(size = 12, color='white',
                                      angle=0),
          strip.background = element_rect(fill = "#525252", color='#525252')) 
  ## make a return obj
  p <- list(fig, df_plot)
  #rename list
  names(p) <- c("plot", " df_plot")
  return(p)
  
}

